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Abstract. We present a new determination of the local dark matter phase-space density. 
This result is obtained implementing, in the limit of isotropic velocity distribution and spher- 
ical symmetry, Eddington's inversion formula, which links univocally the dark matter distri- 
bution function to the density profile, and applying, within a Bayesian framework, a Markov 
Chain Monte Carlo algorithm to sample mass models for the Milky Way against a broad 
and variegated sample of dynamical constraints. We consider three possible choices for the 
dark matter density profile, namely the Einasto, NEW and Burkert profiles, finding that the 
velocity dispersion, which characterizes the width in the distribution, tends to be larger for 
the Burkert case, while the escape velocity depends very weakly on the profile, with the mean 
value we obtain being in very good agreement with estimates from stellar kinematics. The 
derived dark matter phase-space densities differ significantly-most dramatically in the high 
velocity tails-from the model usually taken as a reference in dark matter detection studies, a 
Maxwell-Boltzmann distribution with velocity dispersion fixed in terms of the local circular 
velocity and with a sharp truncation at a given value of the escape velocity. We discuss the 
impact of astrophysical uncertainties on dark matter scattering rates and direct detection 
exclusion limits, considering a few sample cases and showing that the most sensitive ones are 
those for light dark matter particles and for particles scattering inelastically. As a general 
trend, regardless of the assumed profile, when adopting a self-consistent phase-space density, 
we find that rates are larger, and hence exclusion limits stronger, than with the standard 
Maxwell-Boltzmann approximation. Tools for applying our result on the local dark matter 
phase-space density to other dark matter candidates or experimental setups are provided. 

Keywords: dark matter theory, dark matter experiments, rotation curves of galaxies, dark 
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I Introduction 

The direct detection technique has played a major role in the quest for the identification of 
the dark matter (DM) component of the Universe. The goal is measure the recoil energy 
due to the scatterings on a target material of the particles forming the Milky Way DM halo. 
The prime focus is on weakly interacting massive particles ( WIMPs) , non-relativistic thermal 
relics from the Universe with mass in the range, say, between few GeV and few TeV, which 
would scatter elastically on a nucleus, depositing energies of the order of few (tens of) keV (for 
a recent review on the DM problem and direct detection see, e.g., [1]). In recent years several 
experiments have reached the sensitivity to start probing the level of scattering cross sections 
expected for WIMPs. In particular, three collaborations have published results compatible 
with a positive signal: DAM A and DAMA/LIBRA [2] detected an annual modulation in the 
total event rate consistent with the effect expected from WIMPs scatterings because of the 
orbit of the Earth around the Sun; CoGeNT has recently confirmed [3] the detection of a 
low-energy exponential tail in their count rate consistent with the shape predicted for the 
signal from a WIMP slightly lighter than 10 GeV, as already found previously [4], showing 
in addition a 2.8 a indication in favor of an annual modulation effect; finally, CRESST- 

II [5] has just reported the detection of a number events in the acceptance region of low 
energy nuclear recoils, rejecting the hypothesis that they are due to background only, and 
showing the indication for a component of mass either around 25 or 11 GeV. On the other 
hand, other experiments such as CDMS [6] and XenonlOO [7] have not found any evidence 
for a DM signal and produced exclusion plots in the plane WIMP mass versus scattering 
cross section which seem to disfavor the region of the parameter space preferred by DAMA, 
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CoGeNT or CRESST-II (in this respect, see also the recent reanalysis of the data taken at 
an early stage by CDMS at a shallow site [8], and the result of XenonlO [9]). 

The interpretation of a count rate or an annual modulation effect in a given experiment 
in terms of properties of DM particles (and hence obviously the comparison of results ob- 
tained by different collaborations, possibly, with different techniques and/or different target 
materials) depend on a number of assumptions implemented in the analysis of the data. 
Besides issues regarding understanding the target material, such as nuclear form factors or 
whether channeling occurs, and the performance of the detector, such as determining the 
energy threshold, the quenching factors and background rejection/contamination, there are 
uncertainties related to the DM signal itself. First of all there is an uncertainty in the nor- 
malization of the incident DM particle flux, which scales with the local halo density, often 
quoted to be unknown within a factor of 2 or so. Large uncertainties also affect the energy 
spectrum of the DM particles in the detector frame, in turn connected to their velocity distri- 
bution in the Galactic frame and to the proper motion of the Sun/Earth system. Concerning 
these, the vast majority of the analyses adopt a standard paradigm in which the velocity 
distribution is assumed to be Maxwell-Boltzmann with velocity dispersion scaled up of a 
factor of Y^3/2 with respect to the Sun circular velocity (usually taken as the lAU standard 
value of 220 km s~^) and truncated to the value assumed for the escape velocity. 

The Maxwellian distribution is the configuration maximizing the entropy for a self- 
gravitating collisionless system (complete wash out of the initial conditions after gravitational 
collapse) and is associated to the spherical isothermal sphere density profile, which declines 
as at large radii and hence supports a flat rotation curve. It is a well-motivated form 
but unlikely a fair description of the Milky Way DM halo. In particular cosmological N-body 
simulations find that DM halos have density profiles falling more rapidly at large radii, as r~^, 
and velocity distributions showing significant departures from the Maxwell-Boltzmann shape, 
see, e.g., the results from the high-resolution simulations Via Lactea [10] and Aquarius [11]. 
A few analyses have discussed the impact on direct detection results when using DM phase- 
space distribution function as directly read out from the numerical simulations or from other 
specialized approaches devised to describe in detail the fine-grained structure of the DM 
velocity distribution, see, e.g., [12-21]. The main shortcoming of these approaches is that, 
since it is highly non-trivial to include baryons and baryonic feedback in the simulations, 
these analyses treat the Galaxy as if made of DM only (in some cases, normalizing the 
circular velocity obtained in the model for the DM particles to the locally measured circular 
velocity). In reality, the stellar and gas components dominate the potential well in the inner 
Galaxy and out to, at least, our Galactocentric distance, with, e.g., the DM contribution 
to local surface mass density standing essentially within the error bar from star population 
counts [22]. 

Taking an observationally driven point of view, in Ref. [23] (hereafter Paper I) we recon- 
sidered the problem of constructing mass models for the Milky Way, and showed that, when 
considering an updated set of dynamical tracers and an efficient way of scanning a multidi- 
mensional parameter space, a high precision determination of the local DM halo density can 
be obtained, namely around 0.39 GeV cm~^ with a l-cr error bar of about 7%. This result 
relies partially on the assumption of spherical symmetry we imposed on the DM component 
(and on the same prescription to compute the gravitational potential described in section 
2 of the present paper), and would be revised if one argues, for instance, for some level of 
flattening of the DM halo, see, e.g., [24-27]. At the same time, it is a rather solid result 
since there is no tension between the simplified picture we proposed and any of the very var- 



- 2 - 



iegated sets of complementary observables implemented in the analysis (this, together with 
the effectiveness of Bayesian inference applied to our mass model for the Galaxy, is reflected 
in the tiny error we found) . 

In case one assumes that the DM distribution is also isotropic, there is a one-to-one 
correspondence between the spherically symmetric density profile and the underlying distri- 
bution function, with the latter that can be computed from the former through Eddington's 
formula [28]. Although this requires some heavy numerical integrals, it is nowadays possible 
to perform this inversion on very large samples of trial cases. Relying on the same mass 
models introduced in Paper I and an analogous Bayesian approach with a Markov Chain 
Monte Carlo scan of the underlying parameter space, we study here the family of isotropic 
distribution functions one can associate to the spherical DM density profiles. Following this 
approach, we can address here for the first time the theoretical uncertainty on the direct 
detection signal within a framework in which the value of the local halo density, the shape 
of the distribution function and its truncation at the escape velocity, and the circular ve- 
locity of the Sun are taken self-consistently and in agreement with the available dynamical 
constraints. Such method is way more powerful than deriving an overall theoretical error 
assuming that the local halo density, the velocity dispersion, the escape velocity and the Sun 
circular velocity have given uncertainties to be propagated as uncorrelated errors, as actually 
done in most analyses in the literature. 

Results presented in this paper are valid in the limit of spherical symmetry^ and isotropy, 
indeed rather strong assumptions, and should be regarded as a first step towards a study 
allowing for, at least, axisymmetric configurations (the case with axisymmetric models can 
be in principle treated in a specular way, but it is computationally much more demanding; 
such case is subject of ongoing work). On the other hand, the inner regions of a galaxy, such 
as at own position within the Milky Way DM halo, are those for which the simulations find 
weaker evidence for departure from spherical symmetry, and, favored also by the presence of 
large amount of baryonic matter, those with the largest chance for gravitational relaxation 
of the collisionless DM system, and hence where the distribution function is expected to be 
close to isotropic. 

The paper is organized as follows: In Section 2 we introduce the Galactic model and 
the procedure adopted to compute the local distribution function. In Section 3 we review 
the computation of direct detection rates, emphasizing the connection to the distribution 
function, and detail the method implemented to compute the exclusion limit in one sample 
case, the one from the recent data release from the Xenon Collaboration. Section 4 illustrates 
the statistical method implemented, while Section 5 contains our results. Section 6 concludes. 

2 Galactic model and DM phase-space distribution function 

Assuming that the distribution of DM particles in the Galaxy is spherically symmetric and 
isotropic, and in the limit of spherical symmetry for the underlying gravitational potential 
for the Galaxy $(r), Eddington's formula [28] gives an one-to-one correspondence between 
the DM halo density profile Phir) and its phase-space distribution function 



^For the DM profile and for the gravitational potential, which has been "symmetrized" according to the 
prescription presented in section 2. 




(2.1) 



-3- 



where we defined the relative potential ^(r) = — <I>(r) + $(r = Rvir), with the virial radius 
Rmr being the radius at which the DM halo as a virialized object is truncated. Fh depends 
on particle velocity and position in the Galaxy only through the value of the relative energy 
£ = —E + $(r = Rvir) = — -Ekin + "^(r), with E and -Ekim respectively, the total and kinetic 
energy. From the numerical point of view, it is actually simpler to implement Eq. (2.1) by 
changing the integration variable from ^ to the radius of the spherical system r; to compute 
Ffi at the local Galactocentric distance, it is sufficient to specify the radial dependence of ph 
and $ from the local position in the Galaxy out to Rmr- We will then make an ansatz for a 
mass model for the Galaxy, including a parameterization for the contributions of the stellar 
bulge/bar and stellar disc, compare it against available dynamical constraints and implement 
a prescription to extrapolate $(r). 

The Galactic mass model we adopt is the same introduced in Paper I; we summarize it 
briefly here. The DM halo component takes the generic form: 

Ph{r)=p'f{r/ah), (2.2) 

where f{x) is the function defining the shape of the DM density profile. We will consider three 
separate choices, namely: the profile originally proposed by Navarro, Prenk and White [29] as 
the universal profile describing DM halos in numerical N-body simulations, and extensively 
used in the literature: ^ 

fNFw{x) = — — — , (2.3) 

the Einasto profile [30, 31], favored by the latest simulations: 

2 



fsix) = exp 
and, finally, the Burkert profile [32]: 



(x"^ - 1) 
as 



(2.4) 



^^(^)=(l + x)(l + x2)' (2.5) 

which has core, possibly reflecting the wash out of the central DM enhancement induced 
by the baryon infall. In Eq. (2.2), p' and fix, respectively, a reference normalization and 
length scale; most often these are given in terms of the virial mass Myir and the concentration 
parameter Cmr, inverting the relations: 



M,ir = —KnrPO Rt„ = " / r'ph{r) (2.6) 

3 i2DM Jo 

Cvir = Rvir/r-2 (2.7) 

where, in the first equation, which defines also the virial radius Rmr, the virial overdensity 
Amr is computed according to Ref. [33] , while po is the mean background density today. Qbm 
and r^b, the dark matter and baryon energy densities in units of the critical density, enter in 
this expression since we are assuming that only a fraction equal to ilDM/(f^DM + ^h) of the 
total virial mass consists of dark matter; their values are fixed according to the mean values 
from the fit of the 7- year WMAP data [34]. Finally, in the equation for Cyir, r_2 is the radius 
at which the effective logarithmic slope of the profile is —2. 
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For what regards the luminous mass components, the stehar disk mass density profile 
is assumed to take the form: 



Pd{R,z) 



2zd 



sech 



Zd 



with R < Rdn 



(2. 



where is the central disk surface density, R^ and are the length scales in the cylindrical 
coordinates R and z, while Rdm is the truncation radius. Among these, only and Rd will 
be treated as free parameters. Rdm will be assumed to scale with the local Galactocentric 
distance i?o according to the prescription Rdm. = 12 [1 + 0.07(i?o — 8 kpc)] kpc [35]; we will 
also fix the vertical height scale to the best fit value suggested in Ref. [35], Zd = 0.340 kpc, 
since the dynamical constraints we consider in our analysis are insensitive to a slight variation 
around this value. The bulge/bar region is modeled through the mass density profile: 



Pbb{x,y,z) = pbb 



-1.85 



exp(- 



-Sa) + exp 



(2.9) 



where 



]{x'^ + y'^) + z'^ 



and 



+ 



+ 



(2.10) 



The first term in Eq. (2.9) represents an axisymmetric nucleus [36], while the second term 
describes a triaxial bar. When comparing against dynamical constraints, we will implement 
an axisymmetrized version of Eq. (2.9), and assume Xb — Vb = 0.9 kpc • (8 kpc/i?o)i Zb = 
0.4 kpc- (8 kpc/i2o) and qb = 0.6. As for the stellar disk, we fixed these parameters because of 
the lack of observables to discriminate among these values and small deviations around them. 
Such values are in agreement with the best fit obtained in Ref. [37] assuming Rq = S kpc, 
and then scaled to an arbitrary Rq. Rather than using the two mass normalization scales Y^d 
and pbb, in Paper I the parameter space of our reference Galactic model was scanned varying 
the disc mass Md and the balge/bar mass Mbb] here we re-parameterize these in terms of two 
dimensionless quantities, namely, the fraction of collapsed baryons /b and the ratio between 
the bulge/bar and disk masses E, i.e., respectively: 



/b 
E 



OpM + »b Mbb + Md + + Mh^ 

Ob Myir 

Mbb 
Md ■ 



(2.11) 
(2.12) 



In Eq. (2.11) we also included the sub-leading contributions to the total virial mass associated 
with the atomic (Hi) and the molecular (H2) Galactic gas layers, with profiles as given in 
[38], see also Paper I. Einally, in our mass model, the baryons which do not collapse in the 
disc are distributed with the same profile of the DM component. 

The mass density profiles can be used as source terms in a Poisson equation to compute 
the Galactic gravitational potential. A rigorous procedure would require the solution of 
partial differential equations in cylindrical coordinates, a method which would be numerically 
challenging to apply to an extensive scan of our parameter space; moreover, it would actually 
give more information that those required in Eq. (2.1) which assumes a spherically symmetric 
$. In this respect, since the goal is to use Eddington's formula at our position in the Galaxy, 
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i.e. in the outskirts of the stehar Galactic plane, it is a fairly good approximation to infer the 
local deepness of the gravitational potential well, and profile of ^> at larger radii, by applying 
the following prescription^: First, varying the Galactocentric distance we calculate the 
total mass profile M{R) obtained integrating the sum of all density components we introduced 
within a sphere of radius R centered in the Galactic center. This quantity is then used in a 
Poisson equation in the spherically symmetric limit, whose solution <I>(r) reads: 



= Giv 



R^ R 

vir 



(2.13) 



With this last step, one has all the ingredients to estimate the DM phase-space distri- 
bution function according to Eq. (2.1). In particular, since we wish to discuss DM direct 
detection, we will focus on the local DM velocity distribution: 

/o(v) = 1f,,(£:o(v)), (2.14) 
Po 

having factorized out the local DM density po = Ph{Ro) = J d^vFh{£o{v)), and defined 
<So(v) = — ^|vp + ^(i?o)- One quantity which will be useful to compute to compare with 
previous analyses is the local velocity dispersion: 



ai{Ro) = J dS|v|^/o(v), (2.15) 

while one of the advantages of Eddington's approach is that the local escape velocity Vesc is 
not a free parameter, being instead directly related to the gravitation potential through the 
relation: 

v,sc{Ro) = V^R^ ■ (2.16) 

As already mentioned, the DM velocity profile most often used in DM direct detection 
studies is the Maxwell-Boltzmann distribution. This is obtained from the distribution func- 
tion of the isothermal sphere, taking a form in which radial and velocity dependence can be 
factorized: 



Fis{£) oc exp ( ^ j oc pis{r) exp (--^^ ) ^^i |v| < fcsc(r) (2.17) 



and otherwise. Neglecting the small correction due to the truncation at the escape velocity 
(which is assumed here as an external input), this distribution has a constant velocity dis- 
persion, which, after properly normalizing, is found to be equal to a^. The isothermal sphere 
profile has an asymptotically flat rotation curve, with the circular velocity at large radii tend- 
ing to the value 0oo = -\/2/3 u^. Assuming a flat rotation curve for the Milky Way down to 
our position in the Galaxy, that is Qq = Qoo, one usually constrains the Maxwell-Boltzmann 
velocity distribution imposing: 

a^ = y^QQ. (2.18) 



non trivial check of the rehabihty of the method used in this work consists in comparing the velocity 
distribution in the direction perpendicular to the disk, /z, obtained: 1) Within the approximation used in the 
paper. 2) Exactly, computing numerically the axisymmetric gravitational potential associated with our mass 
model for the Galaxy and then performing an Abel inversion to determine fz. We made this comparison for a 
few representative points in our parameter scan and found only a very mild increase in the velocity dispersion, 
and hence a mild impact on the dark matter detection rate. 
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Taking into account the escape velocity in the normahzation factor, the local DM velocity 
distribution in the isothermal sphere case reads: 

1 / 3 \ '^/^ 

^''''^^'^ = N^rX'^v) "'"''^''^'^'^ \A<v...m. (2.19) 

where A^csc is given by: 

iVesc = erf(z) - 2zexp(-^^)/7r^/^ , with z = ^/Jj2vesc/<yv ■ (2.20) 



3 Direct-detection rates and phase-space densities 

The differential direct detection rate for a DM particle x in & given material (per unit detector 
mass) is given by the convolution of the incident flux of the x particles and the differential 
cross section for their scattering off the target nucleus da/dQ [39, 40]: 

- I d'u\u\f^{vi,t) — , (3.1) 



dQ M;vMx^|u|>.„,,, dQ 

where Q is the energy deposited in the detector, Mn the nucleus mass and ■ the x 
phase-space density at the detector and in its reference frame. The lower limit in the integral 
is the minimum velocity required for a x particle to deposit the energy Q. Its expression 
depends from the kinematics of the scattering: In the case of an elastic scattering, such as 
for WIMP-nucleus interactions, the minimum velocity is given by 



QMn 

with Mr = M-^Mn /{M^+Mj\f) the x-nucleus reduced mass, being the mass of the particle 
X- If the scattering is instead inelastic, such as for inelastic DM, the minimum velocity also 
depends from the mass split 5 between the DM particle and the final state; in this case one 
finds 

In this paper we will mainly focus on spin-independent WIMP interactions. Rewriting 
accordingly the differential scattering cross section in terms of the WIMP-nucleon scattering 
cross section at zero momentum transfer a-^n and the nucleus form factor 

where M„ = M^nip/ [M^-\-mp) is the WIMP-nucleon reduced mass and nip the proton mass, 
the differential rate takes the form: 

^ = 2^^'-^-(^) / P2rMhll , (3.5) 

where A is the mass number of the target nucleus. In the prefactor there are quantities 
depending on the particle physics model, while the integrand depends on astrophysical quan- 
tities only, though the limit of integration contains again a dependence on and . 
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The local x phase-space density is computed assuming that the x particles account for 
the entire DM component and is directly associated to Fh{£): 

p^- f^{M,t) = Fh[£^[^)] , (3.6) 

with ^r{t) = Vobs(i) + u. In the last equation Vobs takes into account the motion of the 
detector frame compared to the Galactic frame in which F/^ is computed. One needs to take 
into account the motion of the Local Standard of Rest (LSR) , namely the rotation of the sun 
around the Galactic center vlsr, the motion of the Earth around the Sun and the Sun's 
peculiar velocity V0^pec: 

Vobs(i) = VLSR + V0,pec + V®(t) (3.7) 

with, in Galactic coordinates (where x is the direction to the Galactic Center, y the direction 
of disk rotation, and z the North Galactic Pole) and using (partially) the notation of Ref. [41]: 

VLSR = (0,60,0), 

V0,pec = (C/0, Vq.Wq) ~ (10, 5.2, 7.2) km/s . (3.8) 

In the computations the motion of the Earth around the Sun has been modeled as in the 
DarkSUSY package [42]. Neglecting here for simplicity the ellipticity of the Earth's orbit, one 
finds [43] 

Ve(t) = F®[ficosa;(t-ti) + e2sinw(i-ti)] , (3.9) 

where = 29.8 km/s and ti = 0.218 is the fraction of the year before the Spring equinox, 
while ii and £2 are the directions of the Earth's velocity at times ti and ti + 0.25 years. 

From the previous equations, one can see that astrophysical uncertainties on the event 
rate computation are all encoded in one function, which we will denote by ^^(u, t) and it is 
defined as follows 

5^(u,t) = |u| y d^Fh [£:o(v)] . (3.10) 

This function is the key object in our analysis. For convenience, we will refer to it as to 
the DM phase-space density function, though it is actually an angular integral of the true 
phase-space density boosted in the detector rest frame and multiplied by |u|. 

3.1 Differential rate and exclusion limit for the XenonlOO experiment 

In a real experiment Eq. (3.5) has to be modified in order to account for experimental limita- 
tions related to the efficiency of the detector and its finite energy resolution. In discussing the 
impact of astrophysical uncertainties on the differential event rate, we will consider a sam- 
ple toy-model case, focussing on specifications which would apply, as a first approximation, 
to the case of the Xenon detector, taking a constant efficiency e = 0.3 and including finite 
energy resolution effects by assuming a Gaussian probability density S^{E,Q) that an event 
depositing an amount of energy Q in the detector is detected with an associated energy in 
the interval [E, E -\- dE] (the same assumptions are implemented for the Xenon experiment, 
e.g. see [44] and references therein). When plotting below observed differential event rates 
in an actual experiment, we will refer hence to the quantity:^ 

poo JO 

^ = ej^ dQaE,Q)^{Q), (3.11) 

^ Since we will not be discussing annual or daily modulation effects, we focus on time averaged differential 
event rates. 
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where 

1 _('^_r)^2 



2vr4,(g) 



with the energy dependent dispersion estimated for the Xenon detector: 



CTXc(Q) = VQ/KeV(0.579KeV) + 0.021Q . (3.13) 

We win also discuss how astrophysical uncertainties affect the estimation of the ex- 
clusion limits in the scattering cross section - particle mass plane for a direct detection 
experiment, focusing, as for the differential event rate, on the sample case applicable to the 
Xenon detector. To do that, however, we need to treat finite energy resolution effects in 
a more refined way than what just sketched above. We will follow the approach discussed 
in [45], and take into account that the quantity measured by the Xenon experiment is the 
number of photoelectrons associated with an event, and not the corresponding recoil energy. 
To translate such a signal in to a recoil energy, one has to use the conversion formula [46] 

v{Q) = 2,.QxQCeff, (3.14) 

which gives the number of expected photoelectrons (PE) as a function of the recoil energy. 
There has been some discussion in the literature on which is the correct form to assume for 
the relative scintillation efficiency of the detector >Ce//. In our analysis we will refer to its 
most recent determination [46] providing an estimate for such a quantity at energies as low 
as 3 KeV, while at lower energies we follow [7] assuming that >Ce// goes logarithmically to 
zero. The differential event rate is then given by the following expression 



n=l 



1 



2na 



IT 

PMT 



(3.15) 



where the convolution with the Poisson distribution of mean v{Q) accounts for the fact that 
a given recoil energy has a probability v{Q)'^ / n\ e~'^^^^ of generating n observable PE. The 
Gaussian weighted sum, instead, describes the finite resolution in the observed number of PE: 
the probability of observing Si PE when the actual number is n is assumed to be Gaussian, 
with dispersion equal to y/na-pwr and fJpMT = 0.5 PE. Following [7], we define as signal 
region the interval in the observed number of photoelectrons between 4 PE and 30 PE. In 
this region the Xenon collaboration has reported three candidate DM events associated with 
recoil energies of 12.1 KeV, 30.2 KeV and 34.6 KeV. To extract an upper bound on the cross 
section ^^^n from the three XenonlOO events, one can calculate, as a function of a^n, the 
probability - denoted by C{a^n) - that a random experiment produces a number of events in 
the signal region equal or larger than four. Then, one can exclude cross sections such that, 
for instance, C{(T-^n) > 0.9. The probability C{a-^n) can be written as the sum of different 
Poisson distributions of mean /_f: 

^Kn)=E^^^^"'^'^^"^ (3.16) 
k=A 

where fJ-icr^^n) = A*s(<7;^n) + ^B is the total number of events in the signal region, including any 
known background term fis- In the case of the XenonlOO experiment, we include a constant 
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background contribution equal to 1.8 events. The signal contribution, instead, is given 



where r/ = 48 x 100.9 [kg days] is the exposure of the XenonlOO experiment. Introducing the 
function: 



where in the last line we emphasized the dependence of the exclusion limit a from the DM 
particle mass. 

4 Bayesian analysis of the parameter space 

Eddington's inversion formula Eq. (2.1) establishes - under the assumptions discussed - a 
direct relation between the local DM phase-space density and the parameter space of the 
underlying Galactic model. In turn, the local DM phase-space density, along with the coor- 
dinate transformation in Eqs. (3.7) and (3.8), defines the factor in the DM scattering rate 
depending on astrophysical assumptions. To address the astrophysical uncertainties on DM 
direct detection we need therefore to study the variability range of the Galactic parameters. 
The model we summarized in Section 2 has 7 or 8 free parameters: the local Galactocentric 
distance Rq; the DM halo parameters Ci,ir and M^r^ plus the exponent in case of the 
Einasto profile; the parameters defining mass profiles for the baryonic components iZ^, /{, 
and r (recall Eqs. (2.11) and (2.12)); the halo star anisotropy parameter /3 (introduced to 
make predictions for the velocity dispersion of halo stars and then compare it to the data of 
[47], see Paper I for further details). To set constraints on this parameter space, we consider 
a broad and variegated sample of dynamical constraints for the Milky Way, including several 
recent results, the same set of observables implemented in the analysis on the local halo 
density of Paper I (to which we refer for a comprehensive description): the proper motion of 
stars in the outer Galaxy, the radial velocity dispersion of halo stars, terminal velocities for 
the rotation curve in the inner part of the Galaxy, the monitoring of stellar orbits around the 
Galactic Center, the peculiar motion of Sgr^d*, the Oort's constants, the total mean surface 
density within \z\ < l.lkpc, the local disk surface mass density and, finally, the total mass 
inside 50 kpc and 100 kpc. 

The scan of the parameter space to address variability regions of the parameters is 
performed implementing a Bayesian analysis, again specular to the one carried out in Paper I. 
We sample the Galactic model posterior probability density function (pdf) implementing a 
Markov Chain Monte Carlo scan of the underlying parameter space. By construction, the pdf 
reflects the change of our prejudice about the most probable values for the model parameters 
g after seeing the data d. It is related to the Likelihood C{d\^ and the prior distribution 
-7r(g) through Bayes' theorem, which states 



by: 




(3.17) 



W{a^n) = C{a^n)-^-Q, 
we can define the XenonlOO upper limit (t{M^ as:^ 

VF(c7(M^)) = 0, 



(3.18) 



(3.19) 



P(g|d) 



/:(^|g)^(g) 

p{d) 



(4.1) 



''This condition defines the frequentist 90% confidence level exclusion limit. In the context of the present 
Bayesian analysis, this relation has to be regarded simply as our definition of the exclusion limit as a function 
of the DM particle mass and Galactic model parameters. 



where p(d) is the Bayesian evidence and in the present discussion simply plays the role of a 
normalization constant. For a detailed description of our choice of the Likelihood and of the 
prior distribution we refer the reader to Paper I. In the present work, we generated 16 chains 
for each DM profile, accumulating (after burn in) approximately 1.5 x 10^ accepted points 
in the parameter space for every considered profile. 

Rather than to the posterior pdf of the Galactic model parameters, in this work we are 
interested in the marginal one or two-dimensional pdf of specific functions of the Galactic 
model parameters. We therefore used in our analysis the fact that the marginal posterior pdf 
for any quantity / depending on the model parameters - denoted by p{f,g\d) - is related to 
the posterior pdi p{g\d) by the relation 

p{f,g\d)=6{f{s)-f)p{g\d), (4.2) 

which follows from the definition of conditional pdf. Expectation values and variances of any 
function / of the model parameters with respect to p{g\d) can be then calculated as follows 

(/) = / /(g) P{s\d) ; aj = if) -iff. (4.3) 

Eqs. (4.2) and (4.3) apply to functions of the Galactic model parameters only. The 
phase-space density, the direct detection differential event rate and the XenonlOO exclusion 
limit are instead all functions of both the Galactic model parameters and a further continuous 
variable. Respectively: the detector rest frame velocity u ^, the observed recoil energy E 
and the DM particle mass M^. To keep the present discussion as general as possible, we will 
denote during this section such an additional continuous variable by x. 

The problem is now that of determining the marginal posterior pdf of objects ( - 
functions of the Galactic model parameters - with an explicit dependence from x, i.e. 
C : X — > Ci^)- III our analysis we compute such pdf's discretizing the variable x (see 
also [48] for a different approach to this problem). We therefore use Eq. (4.2) to calculate 
the marginal posterior pdf for a set of N quantities Q = Ci^i)^ with i = 1, - ■ ■ N, where now 
each Q is a function of the Galactic model parameters only. If the steps of such discretization 
{xi ■ ■ ■ xn} are small enough, this procedure leads to a good determination of the function 
C{x). As a result, to any given point x^ in the range of variability of x, we can associate 
instead of a single number, a marginal posterior pdf, namely the pdf related to through 
Eq. (4.2). 

Finally, given a marginal posterior pdf for the quantities Cij one can calculate from 
Eq. (4.3) the means {Q) and construct the required 7% credibility intervals, namely the 
shortest intervals including the 7% of the total probability. In our analysis we considered 7% 
equal to either 68% or 95%. 

5 Results 

The main result of the present analysis is a Bayesian determination of the local DM phase- 
space density based on the Eddington's inversion formula and a broad sample of dynamical 
constraints for the Milky Way. As already mentioned, this quantity incorporates all the 
astrophysical uncertainties which affect the expected event rate at a DM direct detection 

^Indeed, the phase-space density is also a function of time (see Eq. (3.10)). In all figures shown in this 
section we consider a time average of the the function g^. 
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NFW 



Parameters 


mean 


a 


lower 68% 


upper 68% 


lower 95% 


upper 95% 


ft 


0.25 


0.039 


0.21 


0.29 


0.17 


0.33 


eo [km/s] 


246.3 


8.7 


237.6 


254.8 


228.9 


263.6 


Vcsc [km/s] 


550.7 


11.2 


539.7 


561.7 


528.5 


573.3 


ay [km/s] 


287.0 


5.2 


281.7 


292.2 


276.7 


297.2 


Burkert 


Parameters 


mean 


a 


lower 68% 


upper 68% 


lower 95% 


upper 95% 


h 


0.36 


0.049 


0.31 


0.41 


0.27 


0.46 


eo [km/s] 


244.8 


9.1 


235.7 


253.5 


226.0 


261.5 


t;csc [km/s] 


555.4 


10.5 


545.0 


566.0 


534.9 


575.8 


(T„ [km/s] 


293.6 


5.77 


287.8 


299.3 


282.4 


304.9 


Einasto 


Parameters 


mean 


a 


lower 68% 


upper 68% 


lower 95% 


upper 95% 


ft 


0.23 


0.10 


0.13 


0.34 


0.11 


0.46 


eo [km/s] 


245.8 


7.7 


238.4 


253.1 


229.4 


260.6 


liosc [km/s] 


555.0 


17.6 


536.7 


574.0 


523.5 


589.2 


ay [km/s] 


285.1 


570.3 


279.5 


290.8 


274.0 


296.4 



Table 1. Means, standard deviations and confidence intervals for tlie circular velocity Qq, the escape 
velocity Uesc and the dark matter velocity dispersion ay. We also show our results for the parameter 
fb, which affects the shape of the local gravitational potential and therefore the local dark matter 
phase-space density. Other Galactic model parameters have been already discussed in Paper I. 

experiment. We therefore first focus on this quantity; then, in a few illustrative examples, 
we also show how astrophysical uncertainties translates in to error bars for the predicted 
differential event rate and exclusion limit in the scattering cross section - WIMP mass plane. 
Concerning the exclusion limit, we will concentrate on the last XenonlOO data release, how- 
ever, our results for the phase-space density can be readily applied to any other case. 

5.1 Phase-space density 

To better understand our results for the DM phase-space density, it is instructive to analyze 
first a few related local quantities, namely the circular velocity Qq, the escape velocity v^sc 
and, finally, the DM velocity dispersion a^,. In Table 1 we show means and credibility intervals 
for these quantities, while their marginal posterior pdf are shown in Fig. 1, where each 
quantity has been plotted in three cases, corresponding to the three DM profiles considered in 
this work. As we explained in detail in the previous section, such pdf 's have been determined 
applying Eq. (4.2) to the case of the posterior pdf of the Galactic model parameters sampled 
through our MCMC scan. 

The local circular velocity (Fig. 1, left panel) has been very well reconstructed in all the 
three cases. Indeed the three means of Go are all close to the corresponding experimental 
value which has been implemented in the Likelihood function, i.e. Qq = (245 it 10.4) km/s 
[49]. The minor differences between the three cases are related to the fact that different 
profiles have slightly different mean values for Rq and also different total mass within Rq. 

For what concern the escape velocity (Fig. 1, central panel), we find a mean value 
close to 555 km/s for the Einasto and the Burkert profiles, and to 550 km/s in the NFW 
case. Interestingly enough, although we did not consider any direct constraint on the escape 
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Figure 1. Left panel: marginal posterior pdf for the local circular velocity. Central panel: marginal 
posterior pdf for the local escape velocity. Right panel: marginal posterior pdf for the local velocity 
dispersion. 



velocity, our determination of this quantity is in a very good agreement with the central value 
quoted by the RAVE collaboration [50] (median likelihood of 544 km/s and 90% confidence 
interval in the range 498-608 km/s). We also notice that in the Einasto case the Vgsc pdf has 
a broader shape. This result is related to the larger uncertainties - compared to the other 
cases - in the fraction of collapsed baryons that we find in a scan based on an Einasto profile 
(see Table 1). 

The DM velocity dispersion is connected to the width of the phase-space density as a 
function of the velocity. From Fig. 1 (right panel), we can see that in the case of a Burkert 
profile the mean value of ay is significantly larger that in the Einasto and NFW cases. This 
result can be understood as follows. A spherically symmetric gravitational potential generates 
a velocity dispersion for the local population of DM particles which can be estimated through 
the Jeans equation. One finds 

Phir) Jr dr 

By numerical inspection we find that the first derivative of the gravitational potential is 
larger - in the relevant integration region - for a Burkert profile than in the other cases. This 
is related to the fact that, according to our posterior pdf, this profile is in a better agreement 
with observations when a large fraction of collapsed baryons is considered, as one can see in 
Table 1. A larger value of - with respect to the other profiles - is in turn responsible for 
the larger velocity dispersion which emerges from our marginal posterior pdf^. 

We can now focus on the DM phase-space density. In Fig. 2 (left panel) we show the 
marginal posterior pdf for the DM phase-space density evaluated at a reference velocity in 
the detector rest frame (240 km/s). In the notation of the previous section, we are therefore 
considering the case in which x = u and C(x) = gx{u). The u variable has been discretized 
and the reference velocity corresponds to one of the points of such a N-point discretization, 
i.e. Uk = 240 km/s with /c G {1 • • • , N}. Let us note that the Burkert pdf is systematically 
shifted towards larger velocities. This is a feature which is present at each point of the 
discretization. From this marginal posterior pdf (Fig. 2, left panel) we can determine the 

^ In this work (j„ has been calculated through Eq. (2.15). The Jeans equation is however useful to 
qualitatively understand the results 



(5.1) 



CD 




g^(u=340 km/s) [GcV/cm» (s/km)^ ] u [km/s] u [KeV] 



Figure 2. Left panel: marginal posterior pdf for the dark matter phase-space density at u — 240 
km/s. Central panel: mean dark matter phase-space density. Three profiles and Maxwell-Boltzmann 
approximation with 60 = 220 km/s, Vesc = 544 km/s, = 0.3 GeV/cm^ and ay = \/3/2 9o. Right 
panel: ratio between the same Maxwell-Boltzmann distribution of the central panel and the tails of 
our MCMC distributions. 




Figure 3. Phase-space density "bands" corresponding to 68% and 95% credibility intervals. Left 
panel: NEW profile. Central panel: Burkert profile. Right panel: Einasto profile. 



mean value of g^iuk) and the corresponding 68% and 95% credibility intervals. By proceeding 
in this way for all the N points in the discretization we can plot a phase-space density "band" , 
instead of the usual phase-space density curve. This has been done in Fig. 3 for the three 
profiles under analysis. The red band is associated with the 68% credibility interval, while 
the yellow one corresponds to the 95% credibility interval. This figure makes quantitative 
what we already mentioned in the previous section: given a value of the velocity, say u, we 
can now associate to u a pdf - and not just a number - which accounts for astrophysical 
uncertainties in the determination of the phase-space density at that point. 

In Fig. 2 (central panel) we show the " mean" phase-space densities corresponding to the 
three profiles considered in this paper. In this figure to each point of the velocity discretization 
we associated the mean of the corresponding pdf. It should be now clear that, because of a 
larger velocity dispersion, the Burkert case is characterized by a broader phase-space density. 

We now compare our MCMC predictions with the usual Maxwell-Boltzmann approx- 
imation. As already mentioned, a Maxwell-Boltzmann distribution is identified by two pa- 
rameters: the escape velocity Vcsc and the velocity dispersion ay. Then, as explained in 
section 2, the Maxwell-Boltzmann distribution is normalized imposing = a/3/2 ©q, and 
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Figure 4. Mean MCMC dark matter phase-spaee density vs. Maxwell-Boltzmann approximation 
(see the text for more details). Left panel: NFW profile. Central panel: Burkert profile. Right panel: 
Einasto profile. 

this parameter, together with the sharp cutoff imposed at the assumed value for the escape 
velocity, sets the shape of the high velocity tail of the distribution. In Fig. 2 (right-panel) 
we show the ratio between the distribution function found for an isothermal sphere pro- 
file (Maxwell-Boltzmann distribution times the local density, see Eq. (3.6)) with 0o = 220 
km/s, fesc = 544 km/s, = 0.3 GeV/cm^ and (T^, = y^3/2 9o, and the tails of our MCMC 
phase-space densities. From this figure one can see that different profiles lead to distributions 
behaving differently at high velocities, crossing each other in a non trivial way. Interestingly 
enough, we also notice that, in the high velocity tail, our MCMC mean distributions are 
systematically above the ones obtained in the Maxwell-Boltzmann approximation considered 
in this example (and often adopted in the literature as a benchmark distribution), with a 
mismatch which is significantly larger than the mismatch in the numerical value for the local 
halo density (which is about 0.4 GeV/cm^ for the three profiles considered). 

In Fig. 4 we compare our mean MCMC phase-space densities with the distribution 
functions of different isothermal spheres, one for each profile, constructed setting the required 
Galactic parameters at their mean values according to our MCMC scan (see Table 1). From 
this figure we can appreciate the fact that the phase-space density associated with a Burkert 
profile has a shape which - compared to the other cases - is more difficult to approximate 
with the distribution function of an isothermal sphere. As already pointed out, this result is 
related to the higher velocity dispersion found in the Burkert case. 

Finally, in Fig. 5 we show the two-dimensional marginal posterior pdf in the {Qq^cj^) 
plane and compare the shape of this 2D-distribution with the straight line defined by = 



^3/2 00- As one can see, the correlation pattern between these parameters suggested by the 

data is different from the one described by this linear relation. Indeed, the eigenvalues of the 
covariance matrix for the parameters (Go,o"„) point towards directions which are different 



5.2 Differential rate 

As a first application of our results on the DM phase-space density, we now discuss how 
uncertainties on this quantity affect the expectations for the differential event rate at a 
direct detection experiment. The differential event rate can be calculated using Eq. (3.11) 
and depends on the Galactic model parameters through the phase-space density. It is however 
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Figure 5. 2D marginal posterior pdf in tlic {Qo,<^v) plane. Left panel: NFW profile. Central 
panel: Burkert profile. Right panel: Einasto profile. The blue (dashed) line represents the curve 

also a function of the observed recoil energy E which enters in a convolution integral which 
involves the form factor and in the expression for Wmin- Similarly to what we did for the 
phase-space density, we now discretize the variable E and study a set of quantities defined 
by Cfc = ^ (Ek), where E^ is a generic point in the observed energy discretization. Thus, 
to each energy E^. we can now associate a marginal posterior pdf. The associated credibility 
intervals define the "signal bands" which incorporate the astrophysical uncertainties. 

In this section we present the results that we obtained following this approach for three 
illustrative cases: a 100 GeV WIMP (elastic scattering), a 10 GeV WIMP (elastic scattering) 
and, finally, a 45 GeV WIMP (inelastic scattering). The corresponding results are shown for 
the case of an Einasto profile in Fig. 6. In these plots a given credibility interval is the region 
between two curves of the same type. 

From these plots one can appreciate that for light candidates (Fig. 6, central panel), and 
in general in the case of an inelastic scattering (Fig. 6, right panel), astrophysical uncertainties 
can be extremely important and completely spoil the accuracy of the predicted event rate. 
This phenomenon is related to the fact that such candidates - light or inelastic - probe 
only the high velocity tail of the phase-space density. Therefore, a few per cent change in 
Wesc) for instance, can modify the integral over the phase-space density in Eq. (3.5) by a 
factor > 0(1). Instead, heavier candidates scattering off the target nuclei elastically are less 
sensitive to astrophysical uncertainties because these are sensitive to a much larger portion 
of the phase-space density. Thus, in this case, the uncertainties on the differential rate - 
which are due to the uncertainties on the total area subtended by the phase-space density - 
are much less pronounced and are of the order of ~ 20%. 

In Fig. 7 we compare the differential rates obtained assuming different DM profiles both 
for the elastic and for the inelastic scenario. In these plots all the curves have been normalized 
to the differential rate computed assuming a Maxwell-Boltzmann distribution with Qq = 220 
km/s, Vcsc = 544 km/s and /o^ = 0.3 GeV/cm^. In the case of a 100 GeV WIMP the curves 
shown in the plot differ mostly because the three corresponding profiles have different local 
DM densities. In the light and inelastic WIMP scenario, instead, the distance between the 
three curves also depend on the different high velocity behavior of the phase-space densities 
corresponding to the considered profiles. In all cases the Maxwell-Boltzmann approximation 
leads to lower differential rates. 
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Figure 6. MCMC differential event rate for the Einasto profile. Left panel: means and credibility 
intervals for a 100 GeV WIMP elastically scattering. Central panel: means and credibility intervals 
for a 10 GeV WIMP elastically scattering. Right panel: means and credibility intervals for a 45 GeV 
WIMP inelastically scattering with 5 = 100 KeV. 
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Figure 7. Ratios between the MCMC differential event rates and the ones obtained in the Maxwell- 
Boltzmann approximation with Qq = 220 km/s, Vcsc = 544 km/s and = 0.3 GeV/cm"^. Left panel: 
100 GeV WIMP elastically scattering. Central panel: 10 GeV WIMP elastically scattering. Right 
panel: 45 GeV WIMP inelastically scattering with 6 = 100 KeV. 

5.3 XenonlOO exclusion limit 

Finally, we study how astrophysical uncertainties affect the estimation of the exclusion limit 
in the scattering cross section - mass plane. In this section we focus on the recent XenonlOO 
data, implementing a method based on the Poisson statistics, as explained in section 3. We 
expect, however, similar conclusions, i.e similar credibility intervals, for different experiments 
analyzed with more refined statistical tools. 

In the case of the exclusion limit, we discretize the mass of the DM candidate and 
study the quantities Cfc = ^{M^), where is a point in the discretization of the WIMP 
mass, and a is the upper bound defined in Eq. (3.19). From the marginal posterior pdf's of 
the quantities a{M^) we derived a mean exclusion limit and the corresponding "exclusion 
bands" calculated from the credibility intervals associated with the different a(M^). The 
results are shown in Fig. 8 (left panel), where we focus on the case of an Einasto profile. 
Together with the 95% MCMC credibility intervals we also plot for comparison the exclusion 
limits which one finds in the Maxwell-Boltzmann approximation with Galactic parameters 
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Figure 8. Left panel: The blue (dark) band is the MCMC exclusion limit associated with the 
Einasto profile (95 % C.I.). The green (light) band, instead, is the region between the exclusion limits 
obtained in the Maxwell-Boltzmann approximation with respectively Qq = 229 km/s, v^sc ~ 523 
km/s and = 0.33 GeV/cm^ , and Qq = 261 km/s, Vcsc ~ 589 km/s and — 0.45 GcV/cm^. 
Right panel: Ratios between the MCMC exclusion limits (three profiles) and the one obtained in the 
Maxwell-Boltzmann approximation with Qq = 220 km/s, Wgsc = 544 km/s and p^ = 0.3 GcV/cm^. 

set at their upper and lower 95% C.I. values (see Table 1). In agreement with our results 
for the differential events rates, uncertainties at small masses are rather important while at 
larger masses they are less strong (in this respect the logarithmic scale in the plot might be 
misleading) . 

In the right panel of Fig. 8 we show the ratios between the MCMC exclusion limits 
computed for the Einasto, NFW and Burkert profiles, and the one obtained in the Maxwell- 
Boltzmann approximation with Gq = 220 km/s, Vesc = 544 km/s and = 0.3 GeV/cm'^. As 
for the differential rate, the curves differ because the three profiles are charcterized by different 
local densities and velocity distributions, being the tails of such distributions particularly 
relevant for low mass WIMPs. 

Fig. 8 also shows that in all cases the Maxwell-Boltzmann approximation leads to upper 
bounds for the scattering cross section less strong than in the MCMC computation, in par- 
ticular at low masses where the Maxwell-Boltzmann distribution is artificially truncated at 
an escape velocity Vesc = 544 km/s. 

We also studied how astrophysical uncertainties affect the exclusion limit in the scatter- 
ing cross section - mass plane for the inelastic DM scenario. The results obtained assuming 
an Einasto profile are shown in Fig. 9 (left panel). As for the case of DM elastically scatter- 
ing off the detector nuclei, together with the 95% MCMC credibility intervals we also plot 
the results derived in the Maxwell-Boltzmann approximation with Galactic parameters set 
at their upper and lower 95% C.I. values (see Table 1). As expected from our results on 
the differential rates, in the inelastic DM scenario uncertainties on the exclusion limit are 
considerably larger than in the case in which the DM - nucleus interaction occurs elastically. 
Finally, in the right panel of Fig. 9, we show the ratios between the MCMC inelastic exclu- 
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Figure 9. As in Fig. 8 but with 5=120 KeV. 

sion limits obtained for the Einasto, NFW and Burkert profiles, and the one computed in 
the Maxwell-Boltzmann approximation with Qq = 220 km/s, Vesc = 544 km/s and = 
0.3 GeV/cm^. As already mentioned presenting the analogous results for the elastic DM 
scenario, the curves differ bacause of the different local densities and velocity distributions 
characterizing the three considerd DM profiles. 

6 Discussion and Conclusions 

We have presented a new determination of the local DM phase-space density based on the 
Eddington's inversion formula. In this approach, assuming an isotropic velocity distribution, 
and in the limit of spherically symmetric DM density profile phir) and gravitational potential 
for the Galaxy ^{r), there is an one-to-one correspondence between p^ and the DM phase- 
space density. 

Implementing, within a Bayesian framework, a broad and variegated sample of dynam- 
ical constraints for the Milky Way, and sampling through a MCMC algorithm the posterior 
pdf of the Galactic models parameters, we computed the local DM phase-space density. We 
have then applied this result to study a few illustrative cases of direct detection rates, and 
have discussed the impact on the exclusion limit one can infer from the latest data release 
from the Xenon experiment; the extension of our analysis to other DM models or to other 
experimental configurations is straightforward, and similar pattern are expected. The ad- 
vantage of our approach is that it allows to keep a close inspection on the correlation among 
the astrophysical quantities setting the phase-space density of DM particles (and hence in 
turn the direct detection signals) and propagate uncertainties in a consistent way. 

We studied separately three cases corresponding to three possible choices for the DM 
density profile, namely the NFW, Burkert and Einasto profiles. We have found that, on 
average, the velocity dispersion, which characterize the width in the velocity distribution, 
tends to be larger for the Burkert profile than for the NFW or the Einasto. Studying two- 
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dimensional marginal posterior probability in the plane velocity dispersion versus circular 
velocity, we have verified that for any of the three distributions these two quantities do not 
follow the linear correlation pattern suggested in the isothermal sphere model, and adopted 
in most analyses of direct detection rates. We have computed the distributions of the escape 
velocities in our samples and found mean values which are consistent with the mean value 
obtained from observations of high velocity stars by the RAVE Survey; on the other hand, we 
have found that the shape of the high velocity tails in our distributions is sensibly different 
with respect to what you find assuming a Maxwell-Boltzmann distribution truncated at a 
given value of the escape velocity. 

Concerning the DM direct detection signals, astrophysical uncertainties have the largest 
impact on models for which the scattering rate depends mainly on the tails of the phase- 
space density, namely for light WIMPs or DM candidates inelastically scattering, given that 
in such cases the integral of the phase-space density can be considerably different for different 
Galactic models. We have also compared our result with those obtained with a Maxwell- 
Boltzmann defined by setting circular velocity, escape velocity and local halo density to the 
"standard" reference values adopted in most direct detection studies, respectively, 60 = 220 
km/s, Vesc = 544 km/s and = 0.3 GeV/cm^. Regardless of the assumed DM profile, when 
adopting a self-consistent MCMC DM phase-space density, we have found that the differential 
rates are higher and the exclusion limits stronger than in the standard Maxwell-Boltzmann 
approximation. 

Results in this paper are derived under the hypotheses of spherical symmetry and 
isotropy of the velocity distribution, two rather strong assumptions which may be (slightly) 
inaccurate for describing the DM distribution in the local neighborhood (see the note referring 
to this aspect in section 2). Having demonstrated here the power of the Bayesian approach 
in this analysis, we remark that such approach is general enough to be implemented within 
more general Galactic geometries, where extensions of the Eddington's inversion are required 
in order to relate the phase-space density to the parameter space of the underlying Galactic 
model; a closer inspection on this point will be carried out in future work. 
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A Approximate velocity distributions 

The results presented in this paper rely on the MCMC chains generated during our sampling 
of the posterior pdf and have been derived for specific WIMPs masses, cross sections or 
experimental configurations. To apply our phase-space densities to the study of different 
WIMPs or experimental setups, one can use the following approximation 



In Eqs. (A.l) we are approximating the mean MCMC differential rate with the differential 
rate calculated using the mean MCMC phase-space density. This approximation is very 




(A.l) 
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Figure 10. MCMC signals vs tlic signals computed using the MCMC phase-space densities. Wc plot 
Eqs. (A. 3) for a 10 GeV elastic WIMP (left panel), a 100 GeV elastic WIMP (central panel) and an 
inelastic candidate {6 = 100 KeV) of 45 GeV (right panel). Relative errors are of a few % in regions 
where the signal is significantly different from zero. Results have been presented for the NFW case. 



good for ~100 GeV WIMPs and is accurate at the per cent level for light WIMPs or inelastic 
DM candidates in the range of energies where the corresponding signals are significantly 
different from zero (see Fig. 10). Error bars on direct detection signals can be then estimated 
calculating an upper limit and a lower limit cr^ ^ analogously: 



:^A'fUQ) [ du ((5x> ± = , (A.2) 



where (Tg^ can be computed as in Eqs. (4.3). In Fig. 10 we explicitly checked the validity of 
this approximation by studying how the quantities'' 

\ldR\ _ dEL\ I (±) -(±)| 

A = 2x^i^^^^; A±. = 2x%^^ (A.3) 

depend from the recoil energy. As already mentioned, when the signal is non-negligible, 
the relative errors A and A± are of a few per cent or less. Eqs. (A.l) and (A.2) should 
be therefore useful for a quick estimation of DM signals which depend from the local DM 
phase-space density. Files including time averaged {gx)-, idx) ^ ^Qx ^^"^ ^9x) ^ ^Ug^for the 
three considered DM profiles are provided with the electronic version of the paper. 



B Data files provided witii the paper 

With the electronic version of the paper we provide fifteen data files. They contain different 
tabulations of as a function of the velocity u in the detector rest frame as calculated by 
using Eq. (3.10). There are five files for each profile considered in the paper, corresponding 
to {gx)-, {gx) ^ <^Sx ^^"^ ^^x) ^ '^^9x I'espectively. Each file has a name of the form: 

uDF_halotype_curvetype.dat 

with halotype=nfw,burkert,einasto and curvetype=mean, Isplus, Isminus, 2splus, 2sminus. 
Isplus and Isminus correspond to {gx) + CTg^ and {gx) — (Jg^ respectively, while 2splus and 

^ Indeed, in Fig. 10 we also consider finite resolution eflFects according to Eq. (3.11). 
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Effective priors on the derived quantities 




Figure 11. Effective priors on the local circular velocity, escape velocity, velocity dispersion and the 
value of the function g-^^ evaluated at the representative value of 420 km s~^. 



2sminus are instead associated with {g^ + 2ag^ and {g-^ — 2ag^. The velocity is in km/s 
while g^ is in (GeV/cm'^)(s^/km^) and has been time averaged. The format of the files is the 
one required by the DarkSUSY fmiction dshmuDFnum, which means that the first two lines 
are two empty headers and the data just start from the third line. 



C The role of the priors 

In this appendix we will focus on the impact of the choice of the priors on our conclusions. 
In the present analysis we sampled the posterior pdf using flat priors for the Galactic model 
parameters, however, this does not guarantee that the priors actually imposed on the de- 
rived quantities are also non-informative. To determine the effective priors used for the local 
circular velocity, escape velocity, velocity dispersion and the values of the function g^., fol- 
lowing [51], we sampled the posterior pdf for these quantities with a MCMC scan where no 
constraints were imposed on the underlying Galactic model parameters. The pdf found in 
this way provides an estimation of the effective priors used for the derived quantities. In Fig. 
11 one can appreciate that, though the actual priors on these quantities are not "rigorously 
flat", the obtained results are not dominated by the priors as well, since a high probability 
is a priori assigned to a large range of possible values for the derived quantities. 
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